function [dst_img, varargout] = stretchImageByHist(src_img, bin_num, varargin)
%hs.img.stretchImageByHist 图像直方图拉伸
%
% Syntax:
%   dst_img = hs.img.stretchImageByHist(src_img, bin_num[ ...
%               , 'Nodata', nan ...
%               , 'InRange', [] ...
%               , 'OutRange', [0,1] ...
%               , 'HistLowPerc', [] ...
%               , 'HistHighPerc', [] ...
%               , 'HistHalfPerc', 0 ...
%               , 'Gamma', 1]);
%
%   [__, thre_list] = hs.img.stretchImageByHist(__)
%   [__, thre_list, thre_id_list] = hs.img.stretchImageByHist(__)
%
% Params:
%   - src_img      [required ] [numeric        ] 待拉伸图像
%   - bin_num      [required ] [numeric;scalar ] 直方图bin个数
%   - Nodata       [namevalue] [numeric;scalar ] 无效数据
%   - InRange      [namevalue] [numeric;numel=2] 对输入图像的截断范围
%   - OutRange     [namevalue] [numeric;scalar ] 输出图像的范围, 默认为[0,1]
%   - HistLowPerc  [namevalue] [numeric;numel=2] 直方图低阈值百分比
%   - HistHighPerc [namevalue] [numeric;scalar ] 直方图高阈值百分比
%   - HistHalfPerc [namevalue] [numeric;numel=2] 直方图单边阈值百分比, 默认为0
%   - Gamma        [namevalue] [numeric;numel=2] gamma指数
%   - Show         [namevalue] [logical;scalar ] 显示直方图对比图
%
% Return:
%   - dst_img 拉伸后的图像
%   - thre_list 直方图截断两个阈值
%   - thre_id_list 直方图截断两个阈值对应的索引
%
% Author: 
%   wreng, 2024年1月11日
%
% Update: 
%   wreng, 2024年1月11日
%
% Matlab Version: R2023a

%% 参数检查
arguments
    src_img (:,:) {mustBeNumeric}
    bin_num (1,1) {mustBeNumeric}
end
arguments(Repeating)
    varargin
end

narginchk(2, inf);
nargoutchk(0, 3);

%% 解析 name-value 参数
p = inputParser;
% 无效数据, 默认为 nan
addParameter(p, 'Nodata', nan, @(x)(isnumeric(x)));
% 输入截断范围，默认为最小值和最大值
addParameter(p, 'InRange', [], @(x)(isnumeric(x)));
% 输出范围，默认为 0~1
addParameter(p, 'OutRange', [0,1], @(x)(isnumeric(x)));
% 直方图低阈值百分比
addParameter(p, 'HistLowPerc', [], @(x)(isnumeric(x)));
% 直方图高阈值百分比
addParameter(p, 'HistHighPerc', [], @(x)(isnumeric(x)));
% 直方图阈值(单边)百分比，默认为 0
addParameter(p, 'HistHalfPerc', 0, @(x)(isnumeric(x)));
% 指数
addParameter(p, 'Gamma', 1.0, @(x)(isnumeric(x)));
% 显示直方图对比图
addParameter(p, 'Show', false, @(x)(islogical(x)));

parse(p, varargin{:});
nodata = p.Results.Nodata;
in_range = p.Results.InRange;
out_range = p.Results.OutRange;
hist_low_perc = p.Results.HistLowPerc;
hist_high_perc = p.Results.HistHighPerc;
hist_half_perc = p.Results.HistHalfPerc;
gamma = p.Results.Gamma;
is_show = p.Results.Show;

%%
% 过滤无效数据
dst_img = ones(size(src_img))*nodata;
if isnan(nodata)
    valid_mask = ~isnan(src_img);
    src_img= src_img(valid_mask);
else
    valid_mask = (src_img ~= nodata);
    src_img = src_img(valid_mask);
end

% 对输入数据截断
if isempty(in_range)
    in_range = [min(src_img(:)), max(src_img(:))];
end
src_img(src_img < in_range(1)) = in_range(1);
src_img(src_img > in_range(2)) = in_range(2);

% 设置直方图的统计范围
if isempty(hist_low_perc) && isempty(hist_high_perc)
    hist_low_perc = hist_half_perc;
    hist_high_perc = 1 - hist_half_perc;
elseif ~isempty(hist_low_perc) && isempty(hist_high_perc)
    hist_high_perc = 1;
elseif isempty(hist_low_perc) && ~isempty(hist_high_perc)
    hist_low_perc = 0;
    % hist_high_perc = 1 - hist_high_perc;
end
hist_low_perc  = max(0, hist_low_perc);
hist_high_perc = min(1, hist_high_perc);

% 直方图统计, 统计范围默认为 [in_range(1), in_range(2)]
[counts, bin_location] = hsUtils.getImageHist(src_img, bin_num);

% 获取直方图截断范围
cum_counts = cumsum(counts);
hist_low_th = cum_counts(end) * hist_low_perc;
hist_high_th = cum_counts(end) * hist_high_perc;
hist_low_id = find(cum_counts >= hist_low_th, 1, 'first');
hist_high_id = find(cum_counts <= hist_high_th, 1, 'last');
hist_low_val = bin_location(hist_low_id);
hist_high_val = bin_location(hist_high_id);

% 重新进行图像截断
src_img(src_img < hist_low_val) = hist_low_val;
src_img(src_img > hist_high_val) = hist_high_val;
if (hist_low_val == hist_high_val)
    warning("stretchImageByHist warn: hist_low_val==hist_high_val, bin_num(%d) may be too small\n", bin_num);
end

% 归一化
src_img = (src_img - hist_low_val) / (hist_high_val - hist_low_val);

% 指数拉伸
src_img = src_img.^(gamma);

% 缩放到指定范围
dst_img(valid_mask) = src_img * (out_range(2) - out_range(1)) + out_range(1);

if is_show
    figure;
    bar_obj = bar(1:length(counts), counts, 1);
    bar_obj.FaceColor = 'flat';
    
    color_1 = [253 118 63] / 255;
    color_2 = [104 140 200] / 255;

    bar_obj.CData(1:hist_low_id, :) = ones(hist_low_id, 1)*color_1;
    bar_obj.CData((hist_low_id+1):(hist_high_id-1), :) = ones(hist_high_id-hist_low_id-1, 1)*color_2;
    bar_obj.CData(hist_high_id:end, :) = ones(bin_num - hist_high_id + 1, 1)*color_1;
    hold on;
    axis tight;
    max_val = max(counts);
    plot([hist_low_id, hist_low_id], [0, max_val], 'LineStyle', '--', 'LineWidth', 1, 'Color', [0.8500 0.3250 0.0980]);
    plot([hist_high_id, hist_high_id], [0, max_val], 'LineStyle', '--', 'LineWidth', 1, 'Color', [0.8500 0.3250 0.0980]);
end

if nargout > 1
    varargout{1} = [hist_low_val, hist_high_val];
end

if nargout > 2
    varargout{2} = [hist_low_id, hist_high_id];
end


end
    